Radiation Magnetohydrodynamics Simulation of Proto-Stellar Collapse: Two- Component 

Molecular Outflow 

Kengo Tomida^'^, Kohji Tomisaka^'^, Tomoaki Matsumoto^ , Ken Ohsuga^'^, Masahiro N. Machida^, and 

Kazuya Saigo^ 

ABSTRACT 

We perform a three-dimensional nested-grid radiation magneto-hydrodynamics (RMHD) 
simulation with self-gravity to study the early phase of the low-mass stai^ formation process 
from a rotating molecular cloud core to a first adiabatic core just before the second collapse 
begins. Radiation transfer is handled with the flux-limited diffusion approximation, operator- 
spUtting and implicit time-integrator. In the RMHD simulation, the outer region of the first core 
attains a higher entropy and the size of first core is larger than that in the magnetohydrodynamics 
simulations with the barotropic approximation. Bipolar molecular- outflow consisting of two 
components is driven by magnetic Lorentz force via different mechanisms, and shock heating by 
the outflow is observed. Using the RMHD simulation we can predict and interpret the observed 
properties of star-forming clouds, first cores and outflows with millimeter/submillimeter radio 
interferometers, especially the Atacama Large Millimeter/submillimeter AiTay (ALMA). 

Subject headings: stars: formation — ISM: clouds — radiative transfer — magnetohydrody- 
namics 



1. Introduction 

Radiation transfer plays a critical role in star formation and affects the structure of accretion flow and 
the resulting adiabatic cores even in a low-mass regime. However multi-dimensional radiation hydrody- 
namics (RHD) simulations have been rarely performed due to their high computational cost. Therefore, 
the barotropic approximation, which omits radiation transfer and simplifies the thermal evolution of the 
gas, is widely used in multi-dimensional simulations. However, recent advancement in the computer tech- 
nology and development of numerical techniques enable us to incorporate radiation transfer into a multi- 
dimensional magnetohydrodynamics (MHD) simulation within reasonable computational time, using a mo- 
ment method with simplified closure relations. We performed RMHD simulations of proto-stellar collapse 
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with our newly developed numerical code and clarified the difference between RMHD and MHD with the 
barotropic approximation. 

RHD and RMHD simulation codes are being active l y dev eloped and some of them are used to study 
star formation processes. For example, lYorke & Kaisigl (|l995h sol ved radiation hydrodynami cs by using 
the flux limited diffusion (FLD) approximation on 2D nested-grids. IWhitehouse & Batd (|200a) performed 
smoothed particle hydrodynamics simulations with FLD radiation transfer and showed different thermal 
evolution between mode ls with and without ra diation transfer. Orion AMR (Adaptive Mesh Refinement) 
RHD code developed bylKrumholz et all (120071) is used to s tudy the main accretio n phase of low-mass star 
formation (iQffner et al.ll2009l) and high-mass star f ormation (iKrumholz et al. I I2OO9I) . Recently AMR RMHD 
simulations based on RAMSES were reported by ICommercon et all (120091) . who studied the influence of 
radiation on fragmentation. 

Magnetic fields and self-gravity are also the key physical processes in star formation because they 
assume leading roles in the transfer of angular momentum which dominates the global evolution of a cloud. 
Magnetic braki ng and outflow driven by magnetic fie lds efficiently carry the initial angular momentum away 
from the cloud (lTomisakall2002 : Machida et al. Igooeh . If the cloud rotates very fast, it will fragment through 
self-gravitational instability ( Machida et al. 20081) . Since the structure and fate of the core strongly depend 
on its rotation and the efficiency of angular momentum transport, 3D MHD simulations with self-gravity are 
required. 

Molecular outflo ws from young st ellar objects are universally observed with millimeter / submillimeter 
radio instruments (see lArce et al and references therein). In those observations some of outflows show 



comp licated structure which can be interpreted as consisting of multiple components (ISantiago-Garcfa et al. 



2009h . They are strong evidences that magnetic fields play a crucial role in angular momentum transport 
during protostellar collapse. The argument about the driving mechanism of molecular out flow is not yet 



settle d. We propose that the molecular- outflow is driven via magnetic fields in the first core (|Machida et al 



20061) : ho wever, some groups b elieve that the outflow is entrained by a high velocity jet from the central 
protostar (|Velusamy et al.ll2007h . Next generation radio interferometric telescopes like the Atacama Large 
Millimeter/submillimeter Array (ALMA) have the potential to reveal the on-going star formation sites di- 
rectly and contribute to our understanding about overall star formation processes and such important physical 
mechanisms. 

In this letter, we report the results of a low-mass star formation simulation obtained with our newly 
developed 3D nested-grid self-gravitational RMHD code. We performed very deep nested-grid simulation 
of proto-stellar collapse just before the second collapse. In our simulation, a two-component outflow is 
driven from the first core. In this letter, we concentrate on disscussion about these core and outflow including 
comparison between the RMHD simulation and the barotropic approximation while another paper is under 
preparation. The structure of this letter is as follows: In §2 we describe our numerical method and model 
setup. We show numerical results in §3, and §4 is devoted to discussions and conclusions. 
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2. Method and Model 



We solve 3D self-gravitational RMHD equations on nested-grids. Equations and quantities for the 
radiation are defined in the comoving frame of the fluid. We treat only the zeroth order moment equations 
of radiation transfer using the gray (i.e., equations are integrated over frequency assuming a black body 
spectrum of a loc al temperature) FLD approximation proposed by iLevermore & PomraningI (Il98lh with 
Levermord (tl984j) closure relations. We assume ideal MHD with no resistivity. The equation of state (EOS) 
of the gas is also assumed to be ideal, where the gas mainly consists of hydrogen molecules and the adiabatic 
index 7 i s set to 7 /5 throughout the simu lation, which is better to trace the realistic thermal evolution than 
7 = 5/3 (IStamatellos & WhitworthI 120091) . The gas and the dust are assumed to have the common local 
temperature. 

Here we do not describe the implementation of our code in detail but mention some points a bout the 

(l2005bl lal). Our 



Machida et al 



treatment of radiation transfer. For the MHD and self-gravity parts, see 
method for updating radiation and gas energy is similar to ZEUS-MP code ( Hayes et al.ll2006|) . although our 
operator-splitting and treatment of the MHD and radiation parts ai^e different from theirs. While our MHD 
and self-gravitaty parts are advanced explicitly with a second-order scheme, the radiation part is advanced 
implicitly with the first-order backwai^d Euler method for stable time-integration. Newton-Raphson iteration 
is used to solve the resulting non-linear equations. The sparse matrix appearing here is very large and not 
diagonally dominant because of the strong non-linearity of the system; hence, a fast and robust iterative 
solver is required. We use a stabilized biconjugate gradient (BiCGStab) solver with a multi-color modified 
incomplete LU decomposition (ILU(O)) preconditioner, which is stable and efficient. 

We determine the timestep by the Courant-Friedrich-Levy (CFL) criterion derived for the MHD part. 
All the grids have this common timestep and are advanced synchronously although individual (or asyn- 
chronous) timesteps are generally used in AMR simulations. In usual MHD simulations using individual 
timesteps, boundary values are constructed by time-interpolation of the values in the coarser level. However, 
if we use an implicit time-integrator and far lai^ger timestep than that determined by the CFL condition for 
radiation transfer (i.e.. A? » Af/^y = Ax/c), then this time-interpolation is not adequate, at least in principle. 



use 



We use a compiled ta ble of opacity adopted from Ferguson et al.l(|2005l) and lSemenov et al.l(|2003l) . We 



Semenov et al.l (120031) for the Rosseland mean opacity, and we smoothly combi ne the two for the Plan ck 
mean. That is, we use Ferguson et al.l (|2005h in high temperature (T^ lOO OK) and lSemenov et all (|2003l) in 



low temperature because th e Planck mean opacity of lSemenov et al.l (|2003h is quite lower than other opacity 
ta bles in high temperature (IFerguson et al.ll2005h . We adopt the surface formula with the limiter proposed 
by iHowell & Greenoughl (120031) to evaluate radiation energy flux at the cell interface. 



We take a Bonnor-Ebert sphere of lOK with central density pc = l.O x 10~'^gcm~^ for the initial con- 
dition. Initial rotation and magnetic fields are given uniformly along z-axis, w = O.l/f/j ~ 1.5 x 10"^^ sec"' 
and = 1.1/iG. The outermost boundary values are fixed to keep their initial values and no geometric sym- 
metry is assumed. The number of grid points in each level of nested-grids is 64^ . The size of the finer grid 
is half of the coarser grid, and the finer grid is placed around the center of the simulation box self-similarly. 
The simulation starts with 5 levels of nested-grids, and finer grids ai^e generated adaptively to resolve the 
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local Jeans length with 32 meshes not to induce artificial fragmentation (ITruelove et alJll997l) . This is also 
because insufficient resolution causes the gas entropy to be overestimated near the center of the core. At 
the end of the simulation 18 levels of nested-grids are generated and the finest resolution is A;c ~ 0.009 AU. 
We stop our simulation when the gas temperature of the central region reaches 2000K at which temperature 
hydrogen molecule starts to dissociate and the second collapse begins. 



3. Results 



3.1. Overview 



Figure [T] is a typical 3D bird's-eye view of the first adiabatic core and the outflow ~ 500 yr after the 
first core formation. There coexist two components of outflow: well coUimated fast outflow (associated with 
red magnetic field lines in Figure [T]) and slow outflow with a large opening angle (associated with yellow 
magnetic field lines) The fo rmer is driven by magnetic pressure and the latter by magneto-centrifugal force 
(jBlandford & Payndll982l '). iTomisakal (|2002h shows that the magnetic pressure mode typically appears in 
case of weak magnetic fields, while the magneto-centrifugal mechanism appears in case of relatively strong 
magnetic fields. 

Vertical slices of the outflow scale (level L = 10; upper row) and core scale (level L = 12; lower row) are 
shown in Figure |2l Hereafter we discuss only vertical sUces since the process occurs nearly axisymmetric in 
this case. Although the density distribution shows the complicated structure as a result of MHD processes, 
the temperature distribution is almost spherically symmetric. The gas is heated up by the radiation before 
it enters the shock of the first core. This means that the temperature distribution in the outer region is 
dominated by the radiation from the central hot region. Shock heating up to T ~ 30K at the edge of outflow 
is observed. This picture is considerably different from that of the barotropic approximation used in previous 
studies, in which the temperature is determined only by the local gas density. 

The difference between the RMHD simulation and the barotropic approximation is clearly seen in the 
right column (c) in Figure |2] where the ratio of the gas temperature obtained in the RMHD simulation to 
the barotropic temperature given by the local gas density is plotted. In the RMHD simulation, the gas tends 
to attain the temperature typically 2-3 times higher than that in barotropic approximation around lOAU 
from the center, although the two models are very close in the innermost region of the core (<0.1AU). 
The difference is most striking just above the shock at the surface of the first core because of pre-shock 
heating by radiation from the core. Here we use the following polytropic relation to evaluate the barotropic 
temperature: 



r= 10 



Pc 



7-1 



K, 



(1) 



where pc = 2.0 x 10~^^gcm~^ is the critical density and 7 = 7/5 is the adiabatic index of diatomic molecular 
gas. The parameters used h ere are chosen to trace the th ermal evolution track of the central region in a 
spherical RHD simulation by lMasunaga & Inutsukal (l2000h . 
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3.2. Two-Component Outflow 

We visualize the outgoing mass flux and angular momentum flux defined as |pvp| and |pVp|r x vt|| 
where v,- > respectively where Vp(vt) is poloidal (toroidal) velocity at the early driving phase (~ 350 yr 
after the first core formation) of outflow in Figure [3] In the figure we can recognize two components 
of outflow driven by different mechanisms: the outer magneto-centrifugal mode and the inner magnetic 
pressure mode. 

The inner magnetic pressure mode has a relatively high velocity ~ 2 km s~^ and is well collimated. The 
front of this outflow gets hot through shock heating due to its high Mach number. However, the outer, slow 
outflow has the mass and angular momentum flux, one order of magnitude larger, since the gas in the inner 
region has a relatively small angular momentum. This trend is enhanced as the outflow evolves. Accordingly 
we can say that the outer magneto-centrifugal mode is more important for the kinematic structure of the 
accretion flow as a carrier of angular momentum; however, the fast, inner outflow will stand out more for 
observation because it travels a long distance and achieves a high temperature. 



3.3. Structure of the First Core 

We can see a more detailed structure of the collapsing cloud from Figured where we plot the distribu- 
tion of physical quantities at the end of the simulation, just before the second collapse starts (~ 650 yr after 
the first core formation). Panels (a)-(c) show the distribution of physical quantities in the disk midplane and 
along the rotational axis. At this time the fast outflow reaches z ~ 85 AU. The surface of the core is located 
at r ~ 35 AU in the disk midplane and z ~ 4AU along the rotational axis . This height of th e first core is larger 
than that with the barotropic approximation, which is typically ~ lAU (|Saigo et al because the outer 



region of the first core achieves a higher entropy in RMHD simulation than in the barotropic approximation. 
The core with a higher entropy can support more mass even with the same central density; therefore, the 
first core lives longer in RMHD under the same accretion rate, which is mostly determined by the initial 
condition. Along the z-axis there is a high entropy radiative precursor outside the core (z ^ 5AU) seen in (b) 
and (c). The central region within 2AU from the center seems nearly spherically symmetric due to the 
slow rotation because of efficient magnetic braking. 

The structure of the disk and the behavior of angular- momentum transport is visible in panel (d) of 
Figure |4] and Figure |2l The profiles of infalling velocity -v^ and rotational velocity ai^e plotted in panel 
(d) of FigurelH The infalling gas initially decelerates at the weak shock near ~ 35 AU associated with the first 
core. There exist two centrifugal barriers, and the outer one near ~ 20AU con^esponds to the driving region 
of the slow and wide outflow driven via the magneto-centrifugal force. This outflow carries off angular 
momentum of the gas efficiently and the gas falls radially. Nonetheless in this case the initial magnetic 
fields are not so strong, the gas spins up again as it falls and hits the inner centrifugal barrier near ~ 5AU. 
Then, the magnetic field lines ai^e wound up tightly (also due to weak initial magnetic fields) and magnetic 
pressure launches the fast, well collimated outflow. 
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The variety of the gas's thermal behavior is shown in the p-T plane in Figured The distribution 
in this plane indicates a strong dependence on the geometry. We can see again that the gas tend to be 
hotter in the RMHD simulation th an in the barotropic approximation. Although this was already pointed 
out by IWhitehouse & Batd (120061) . our simulation can resolve more detailed features. The jump around 
p ~ 10~'^gcm"^ is caused by the shock at the surface of the first core. This density is very close to the 
critical density in the barotropic approximation. This shock is nearly isothermal since radiation can transfer 
thermal energy produced at the shock to the upstream efficiently. Radiation transfer is not so efficient 
beyond this critical density and this forms the first adiabatic core. When the temperature reaches 1500K 
the evolution track in the p-T plane gets shallow. This is because all the components of the dust evaporate 
and the opacity drops sharply at this temperature, and then radiation transfer becomes effective again. The 
structure of the first core is slightly modified by this effect but not drastically since the gas soon becomes 
optically thick and adiabatic again. We emphasize that radiation transfer is required to handle the realistic 
thermal evolution described here. 



4. Conclusions and Discussions 

We performed numerical simulations of the early phase of the low-mass star formation process just 
before the second collapse starts. Our simulation used a newly developed three-dimensional nested-grid 
self-gravitational FLD RMHD simulation code without assuming artificial thermal evolution. In our case, 
radiation transfer does not seem to change the qualitative scenario of low-mass stai^ formation drastically. 
However, the temperature distribution is significantly changed by introducing radiation transfer. Realistic 
treatment of gas thermodynamics alters some properties and structure of the core quantitatively. The mass 
and size of the first core at a certain central density are larger because of a higher entropy, and the lifetime 
of the first core becomes slightly longer. Furthermore the envelope and the outer region of the first core 
become hotter. We suggest that the observational probability of such very young star forming sites rises 
when compared to previous predictions based on the simulations with the barotropic approximation. 

Barotropic approximation fails to reproduce realistic thermal evolution; needless to say radiation trans- 
fer plays a critical role there even in the low-mass regime. In the barotropic approximation all the gas 
elements trace the thermal evolution track of the gas at the center of the cloud which experiences no shock 
and, therefore, has minimum entropy in the spherically symmetric simulation. In this sense barotropic ap- 
proximation tends to underestimate the gas entropy and temperature. In 3D, the RMHD simulation has more 
striking and complex differences from the barotropic simulation. On the other hand if there exist fast initial 
rotation and ineffective angular momentum transfer, then the entropy of the resulting thin disk-like first core 
may be lower than in the barotropic EOS. Thus radiation transfer will affect the mechanical properties such 
as stability or fragmentation through the thermal property of the gas. Determining a realistic temperature 
distribution is essential for predicting or interpreting the properties acquired in both the molecular emission 
lines and thermal continuum observations, for strong temperature dependences of the emissivity. Chemical 
reactions are also sensitive to temperature, so our results can be applied to a study of chemistry in molecular 
clouds and star forming regions. 
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In our case, unlike previous MHD or RMHD simulations, there coexist two components of outflow 
from the first core. This feature suggests that the angular momentum transfer may occur necessarily not 
in single but in multiple stages, at least under a certain condition of magnetic fields although the situation 
depends on the initial condition and will be changed when magnetic dissipation is incorporated. Nonetheless 
our results imply a rich diversity of the outflow and first core. 

Though on-going star formation like this simulation has never been observed directly, new telescopes 
are expected to reveal them, and ALMA see ms to be the most prom ising in the near future. We propose 



some VeLLOs (Very Low Luminosity Objects: iDunham et al.l (120081) ) are candidates for such obscured star 



forming cores. Compact (~ 100 AU scale) and warm (~ 30 K) molecular outflow can also be a good indi- 
cator of the first core. Our results can be compared with future observation. 
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Fig. 1. — 3D structure of the first core and outflow in level L- 10. The left and bottom panels are density 
profiles and the right panel shows the temperature distribution. The cyan surface is a density isosurface. The 
fast outflowing region (Vj > 0.3 km s~^) is also visualized with volume rendering. Red and yellow lines are 
the magnetic field lines. 
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(a) log (l^s Density ^/ccj) (b) lag. (Gas Temperature [K]) (c) Gas Temperature / BarotrofMC Temperature 



Fig. 2. — Vertical slices of (a) gas density, (b) gas temperature, and (c) gas temperature to the barotropic 
temperature ratio. The upper row is level L = 10 corresponding to a scale of ~ 150AU and the lower row 
is level L = 12 corresponding to a scale of ~ 40AU. The projected velocity of the gas is overplotted with 
arrows. 
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Fig. 3. — Outward mass flux (left) and outward angular momentum flux (right) in level L = 11 at the early 
driving phase of the outflow are graphed. Two components of the outflow are clearly observable. 
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Fig. 4. — Plots of (a) gas density, (b) gas temperature, and (c) gas entropy (K = P/ p^) in the disk midplane 
(r; red solid line) and along the rotational axis (z; blue dashed line), and (d) infalling/rotational (-v,-; red 
solid line / v^; blue dashed line) velocities in the disk midplane at the end of the simulation. 



-13- 




Gas Density [g/cc] 



Fig. 5. — Distribution of thermal properties of the gas at just before the second collapse starts (~ 650 yr 
after the first core formation) in the p-T plane in the disk midplane (r; red crosses) and along the rotational 
axis (z; blue squai^es). The green dashed line is the barotropic EOS used in previous simulations. 



